% T1_mCINE_IR_NoiseEstimation: Estimates parameters and errorbars for T1 recovery in an IR sequence 
% Used signal model: S = A-B*exp(-R1*TI) ] 
% Parameter A = fully relaxed signal, B = inversion efficiency, R1 = 1/T1  (See McKenzie et al. 2006)

% After executing a file browser pops up, where the image to be mapped
% needs to be selected. All images should be in the same folder. The
% selected image should be from the first cardiac cycle.

% a roi where the mapping should take place has to be drawn. The mean
% signal is calculated in that roi, and that is used for the mapping

% There is a variable 'repeatphase', which should be the number of images
% after which the same cardiac phase is reached in the CINE acquisition.

% As an output the estimated values for R1, T1*, A and B are given. If the
% variable txtfile = true, a .txt file is created with T1, STDCRT1,
% Heartrate and Repetition time.

% Created by Henk Smit, EMC, 05-2011 based on the work by Dirk Poot,
% University of Antwerp, 13-8-2007.

clear

% switch for voxel/roi/roimean type of estimation
% WORKS ONLY WITH 2 FOR NOISE ESTIMATION SO FAR, SO NOT USED
% roibased = 2; % 0=voxel, 1=roi all voxels, 2=roi mean
 
% fill after how many images the same cardiac phase is reached
repeatphase = 1;

% which should be the first point to take into account
fitfrom = 1; 

% output values as a txt file
txtfile = true;

[file,chosendir]=uigetfile('*.dcm');
imagenr = str2num(file(1:end-4));
cd(chosendir);
d=dir;
nfol=(length(d)-2)/repeatphase;
cd ..
outputdir = pwd; %save maps in dir of chosen image

%read in the image that is picked by the user
chosenimage = dicomread(dicominfo(fullfile(chosendir,file)));

%get the dicom info
for i=1:nfol   
    %henk orig fitdata(i).folder=char(d((i-1)*repeatphase+2+imagenr).name); 
    fitdata(i).folder=char(d((i-1)*repeatphase+2+1).name); 
    info = dicominfo(fullfile(chosendir,fitdata(i).folder));
    fitdata(i).image = double(dicomread(info));
    fitdata(i).TE=info.EchoTime;
    fitdata(i).TR=info.RepetitionTime+99999; %TR is not correctly in header, should look it up. It's long compared to T1 anyway, so almost neglectible
    fitdata(i).Angle=info.FlipAngle;
    fitdata(i).info=info;
    fitdata(i).TI=info.TriggerTime;
end

%define masks
for i=1:nfol
    if i==1
%      henk orig imshow(fitdata(2).image,[90 570]);
% disp('startin g figure');
%        imshow(double(chosenimage),[])
%        
        figure
        mx = max(chosenimage(:));
        mn = min(chosenimage(:));
        image((double(chosenimage)-double(mn))/double(mx-mn)*256);
        colormap(gray(256));
        axis equal
        set(gca,'xtick',[],'ytick',[]);

% % when ROI is to be drawn here:
% % disp('starting roipoly');
       M=roipoly;
% % disp('done roipoly');
       close;
       
% To import ROI:
%       load('C:\DicomQRPush\mr_physics_t1_mapping_mcine_2086\maskphantom196.mat')

% To save ROI:
%       save('maskphantom1962.mat','M')

    end

    fitdata(i).image=fitdata(i).image.*M;

end

for k=1:nfol
    TI(k)=double(fitdata(k).TI);
    TR(k)=double(fitdata(k).TR);
end

%initialize&allocate
TI=TI';
TR=TR';

yd=mean(double(nonzeros(fitdata(1).image)));

for k=2:nfol
    yd=[yd mean(double(nonzeros(fitdata(k).image)))];
end

yd=yd';
nrvoxels = 1;      

% LS fit for initialization
fiteq2 = ('a-b*exp(-x*c)');
s2 = fitoptions('Method','NonLinearLeastSquares','Lower',[0,0,0],'Upper',[30000,30000,30000],'Startpoint',[400, 800, 0],'Maxiter',500,'Display','off','TolFun',10^-15,'TolX',10^-15);
f2 = fittype(fiteq2,'options',s2);
[valmin,locmin]=min(yd);
lsstart = locmin+1;
yd2 = yd(lsstart:end);
TI2 = TI(lsstart:end);
[c22,gof,output] = fit(TI2,yd2,f2,s2);
A2=c22.a ;
B2=c22.b;
R12=c22.c;
initval = [1000*R12, A2, B2, 30]';

initval = [1, 500, 1000, 30]';
initR1s=0.3:0.3:7.5;
initAs=200:75:1800;
initBs=300:75:3000;
initvalmat=[0.45 300 400 30];

for rr=1:size(initR1s,2)
    for aa=1:size(initAs,2)
        for bb=1:size(initBs,2)
            initvalmat = [initvalmat; initR1s(rr) initAs(aa) initBs(bb) 30];
        end
    end
end
initvalmat=initvalmat';

extraInitialValues = mat2cell(initvalmat,4,ones(1,size(initvalmat,2)));


% Estimation
fun = @(tht) predict_IRT1( tht , TI(fitfrom:end));
tht = fit_MRI( fun, yd(fitfrom:end), initval,'numPDFoptpar', 1,'initialValueSpecifierVect',[1 0 0 1],'extraInitialValues',extraInitialValues); 
[CRLB, I, J] = CramerRaoLowerBound_MRI( tht(1:3,:,:), fun, tht(4,:,:));

ML(1) = tht(2);
ML(2) = tht(3);
ML(3) = tht(1);
ML=ML';
    
A = ML(1);
B = ML(2);
T1 = 1000/ML(3);
R1 = ML(3)/1000;
CRLBsigma = (size(yd,1)-fitfrom)*tht(4)/(size(yd,1)-size(tht,1)-fitfrom+1);
CR = zeros(3,3);
CR(3,3)=CRLB(1); 
CR(1,1)=CRLB(3);
CR(2,2)=CRLB(6);
    
STDCRR1 = sqrt(CR(3,3)/1000);
STDCRT1 = sqrt(CR(3,3))/(1000*R1^2); % c*A = c*std(A)
STDCRA = sqrt(CR(1,1));
STDCRB = sqrt(CR(2,2));


%display plots and results

scatter(TI,yd,'*r')
hold on;
ti=floor(min(TI)):1:round(max(TI));
yy=abs(A-B.*exp(-ti*R1));
yy2=abs(initval(2)-initval(3).*exp(-ti*0.001*initval(1)));
plot(ti,yy,'--b');
xlabel('Inversion Time (ms)')
ylabel('Signal Intensity')
disp(['Model = Abs(A - B*exp(-R1*TI))'])
disp(['Results: value +- standard deviation:'])
disp(['R1 = ' num2str(R1) ' +- ' num2str(sqrt(CR(3,3)))])
disp(['T1* = ' num2str(T1) ' +- ' num2str(STDCRT1)])
disp(['A = ' num2str(A) ' +- ' num2str(STDCRA)])
disp(['B = ' num2str(B) ' +- ' num2str(STDCRB)])
disp(['T1, A and B are saved in ' fullfile(pwd,'mCineEstimationOutput.txt')  ])

%write to txt file if desired
if txtfile == 1;
    outputfile = fullfile(pwd,'mCineEstimationOutput.txt');
    fid = fopen(outputfile, 'a');
    fprintf(fid,'%i %i %i\n',T1, A, B); 
    fclose(fid);
end

% for executable add:
% waitfor(gcf)
% 
% exit;    
